i = 0
nn.results <- list()
size_seq <- seq(10,20,2)
# Loop over years
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Get optimal parameters from cross validation
  # Fit Neural Network
  nn.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("nnet")) %dopar% {   # Loop over CV runs
                           set.seed(52*cv) 
                           shuffle <- sample(1:N,
                                             N,
                                             replace=FALSE)
                           predictions <- matrix(NA,nrow=N,ncol=length(size_seq))
                           sampleSplit = c(0,round((N/5*c(1:5))))
                           
                           for (f in 1:5) {
                             # Remove constant variables from predictors
                             train.RHS <- dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                        rhs]
                             test.RHS <- dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                       rhs]
                             # Remove variables that are constant
                             not.constant <- apply(train.RHS,2,sd)!=0
                             train.RHS <- train.RHS[,not.constant]
                             test.RHS <- test.RHS[,not.constant]
                             # Find PC weights
                             pcs.fit <- prcomp(train.RHS,center = TRUE, scale = TRUE)
                             train.pcs <- predict(pcs.fit,
                                                  train.RHS)[,1:min(nn_count_params$pcNum,
                                                                    ncol(train.RHS))]
                             test.pcs <- predict(pcs.fit,
                                                 test.RHS)[,1:min(nn_count_params$pcNum,
                                                                  ncol(train.RHS))]
                             # Fit Neural Network
                             for (s in 1:length(size_seq)) {
                               # Fit Neural Network
                               cv.fit = nnet(x = train.pcs,
                                             y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                         dv]),
                                             size = size_seq[s], 
                                             decay = nn_count_params$decay,
                                             trace = nn_count_params$trace,
                                             linout = nn_count_params$linout
                               )
                               predictions[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                         s] <- predict(cv.fit,
                                                       newdata = as.matrix(test.pcs))
                             }
                           }
                           predictions[predictions<0] <- 0 
                           cv_results <- c("CV Run"=0,
                                           "size"=0,
                                           "R2"=0)
                           for (s in 1:length(size_seq)) {
                             R2 <- Rdev(dta.past[,dv],
                                        predictions[,s])
                             print(R2)
                             if (R2>cv_results["R2"]) {
                               cv_results <- c("CV Run"=cv,
                                               "Size" = size_seq[s],
                                               "R2" = R2)
                               best.s <- s
                             }
                           }
                           return(cv_results)
                         } 
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  # Fit Neural Networks
  nn.results[[i]]$size <- round(median(cv_results[,"Size"]))
  dta.past.rhs <- dta.past[,rhs]
  
  not.constant <- apply(dta.past.rhs,2,sd)!=0
  dta.past.rhs <- dta.past.rhs[,not.constant]
  pcs.fit <- prcomp(dta.past.rhs,
                    center = TRUE, 
                    scale = TRUE)
  insample.pcs <- predict(pcs.fit,dta.past.rhs)[,1:min(nn_count_params$pcNum,ncol(dta.past.rhs))]
  outsample.pcs <- predict(pcs.fit,dta[dta[,t.var]== year,])[,1:min(nn_count_params$pcNum,ncol(dta.past.rhs))]
  
  set.seed(i)
  nn.results[[i]]$mod <- nnet(x = insample.pcs,
                              y = dta.past[,dv],
                              size = nn.results[[i]]$size,
                              decay = nn_count_params$decay,
                              trace = nn_count_params$trace,
                              linout = nn_count_params$linout
  )
  nn.results[[i]]$fit.oos <- predict(nn.results[[i]]$mod,
                                     newdata = outsample.pcs)
  nn.results[[i]]$fit.oos[nn.results[[i]]$fit.oos<0] <- 0 # Left censor predictions
  nn.results[[i]]$actual.oos <- dta[dta[,t.var] == year,
                                    dv]
  print(year)
}
save(nn.results,
     file=paste(modeldir,"/",
                country,
                "_nn_",
                "count",
                "_",
                rhs.group,fileext,
                ".RData",
                sep = ""))